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We reformulate the coupled set of continuum equations for the renormalized gluon and ghost 
propagators in QCD, such that the multiplicative renormalizability of the solutions is manifest, 
independently of the specific form of full vertices and renormalization constants. In the Lan- 
dau gauge, the equations are free of renormalization constants, and the renormalization point 
dependence enters only through the renormalized coupling and the renormalized propagator 
functions. The structure of the equations enables us to devise novel truncations with solutions 
that are multiplicatively renormalizable and agree with the leading order perturbative results. 
We show that, for infrared power law behaved propagators, the leading infrared behavior of the 
gluon equation is not solely determined by the ghost loop, as concluded in previous studies, but 
that the gluon loop, the three-gluon loop, the four-gluon loop, and even massless quarks also 
contribute to the infrared analysis. In our new Landau gauge truncation, the combination of 
gluon and ghost loop contributions seems to reject infrared power law solutions, but massless 
quark loops illustrate how additional contributions to the gluon vacuum polarization could re- 
instate these solutions. Moreover, a schematic study of the three-gluon and four-gluon loops 
shows that they too need to be considered in more detail before a definite conclusion about the 
existence of infrared power behaved gluon and ghost propagators can be reached. 
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In the Standard Model of the strong, weak and electromagnetic forces, the interactions are quantitatively de- 
scribed by gauge field theories. Quantum chromodynamics is a non-Abelian gauge theory, and the proof of its 
renormalizability [jl] and discovery of ultraviolet asymptotic freedom^ have been milestones in its acceptance as the- 
ory of the strong interaction. For large momenta, perturbation theory seems an appropriate calculational tool, as the 
1} . coupling becomes small. However, for small momenta the coupling grows large and adequate methods have to be used 
'■^ ' to study nonperturbative phenomena like confinement, chiral symmetry breaking and fermion mass generation j3|. One 
such method is the study of the Dyson-Schwinger equations Q, and their phenomenological applications to hadronic 
physics is a subject of growing interest^]. The gluon self- interaction and existence of ghost fields are remarkable 
features of non-Abelian gauge theories, and the study of the infrared behavior of the gluon and ghost propagators, 
■ and of the running coupling has become a major part in the quest to understand color confinement in QCD. 

Early studies of the Dyson-Schwinger equation for the gluon propagator in the Landau gauge concluded that the 
gluon propagator is highly singular in the infrared^, 0, ||, However, these studies neglected any contribution of 
the ghost fields, and moreover, required one to assume cancellations of certain leading terms in the equations. It is 
therefore far from certain that these solutions have the correct QCD infrared behavior. That these solutions are at 
the origin of successful phenomenological applications can be understood by observing that they can generate the 
necessary integration strength in the kernels of the gap and Bethe-Salpeter equations JiT| . 

Recent studies of the coupled set of continuum Dyson-Schwinger equations for the renormalized gluon and ghost 
propagators, initiated by von Smekal, Hauck and Alkofer|ll|, and later confirmed by Atkinson and Bloch[[l2j [l3|], have 
shown that the coupling of ghost and gluon fields plays a crucial role in the generation of a consistent infrared behavior 
of QCD. In the infrared the gluon and ghost propagators are power behaved, and the strong running coupling has an 
infrared fixed point. The leading infrared exponent is determined by requiring the power solutions to be consistent 
with both gluon and ghost equations. Although its precise value depends on the vertex Ansatze and other details of 
the truncation, all the results showed that the ghost propagator is more divergent than its bare counterpart, while the 
gluon propagator is less divergent, and even vanishes as the momentum goes to zero. Furthermore, another common 
feature of these studies is that the leading infrared power in the gluon equation depends only on the ghost loop, and 
not on the gluon loop. One of our aims is to investigate whether this is a characteristic of QCD or merely an artifact 
of the approximations. 

Multiplicative renormalizability is an important feature of gauge field theories, and it is ruined by the approximations 
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to the renormalized vertices and to the vertex renormalization constants introduced in these studies. Furthermore, the 
high momentum behavior of the approximate solutions is generally not in agreement with the results of perturbation 
theory. Herein we aim to develop a truncation scheme, which respects multiplicative renormalizability, and has the 
correct perturbative limit as well. 

Let us first briefly discuss the approximations introduced previously and the defects of their solutions. In Refs. |l^, 
|l3f bare vertices are used and the triple-gluon vertex renormalization constant Z\ is set to its bare value of one. This 
simple approximation allows a straightforward basic analysis of the coupled set of ghost and gluon equations. As has 
been shown in detail in Ref. fl2|| , this approximation violates the multiplicative renormalizability of the solutions, 
and furthermore, neither the anomalous dimension of their logarithmic ultraviolet behavior nor the leading order 
beta coefficient of the running coupling satisfy the results of perturbation theory. The authors were, however, not 
able to conclude if this slight, but apparent, discrepancy with perturbation theory was due to an artifact of their 
approximation or if it could be a genuine feature of full QCD. In Ref. [ pT| more care was given to the construction of 
the vertices, requiring the Ward-Takahashi identities to be satisfied to some extent. However, to force the solutions 
to have the correct perturbative behavior, the authors introduced a manipulation of the renormalization constant Z\, 
which again violates the multiplicative renormalizability of the solutions. 

In Sections (III- VIII) of Ref. Jl2[ , an intriguing truncation was studied: based on symmetry ideas between the 
ghost and gluon equations only the ghost loop is retained in the gluon equation. The paper showed in detail that 
the solutions of this truncated set of equations are multiplicatively renormalizable and that their high momentum 
behavior satisfies the perturbative results of this truncated theory. Although this truncation is not consistent with 
QCD, as the gluon loop was discarded, its solutions exhibit these important properties expected from solutions of gauge 
field theories, thereby providing insight into how multiplicative renormalizability of the nonperturbative solutions is 
achieved, how the high momentum behavior of these solutions obey the leading order perturbative results, and even 
how both features are entangled. Section (IX) of Ref. |l^| then went on to show how a naive treatment of the gluon 
loop destroys both multiplicative renormalizability and agreement with perturbation theory. 

In addition to the above-mentioned approximations, Refs. jO, [l2| also introduce angular approximations to make the 
analytical and numerical calculations more tractable. These angular approximations are satisfying for slowly varying 
functions , like propagator functions in their perturbative regime. However, Ref. showed that these angular 
approximations can introduce a significant error for regions where the functions are power behaved. Although this 
means that eventually the analysis of the asymptotic infrared power behavior will require the angular integrations to 
be performed exactly, the qualitative conclusions about multiplicative renormalizability and high momentum behavior 
of the solutions are, however, not affected by the angular approximation. 

In this paper we reformulate the coupled set of continuum equations for the renormalized gluon and ghost propaga- 
tors, such that the Landau gauge equations are free of renormalization constants, and depend on the renormalization 
point only through the renormalized coupling and the renormalized propagator functions. In analogy to the ghost-loop 
only truncation of Ref. p"2[ , we see how the simple structure of the equations makes the multiplicative renormalizability 
transparent, i.e. not hidden in renormalization constants or renormalized vertices. We then introduce a novel trunca- 
tion scheme to these formal equations, which respects both the multiplicative renormalizability and the perturbative 
limit of the solutions. We show that, for power law solutions, the leading infrared behavior of the gluon equation no 
longer solely depends on the ghost loop, but rather depends on the ghost loop, the gluon loop, the three-gluon loop, 
the four-gluon loop and even on massless quark loops. The consequences for the existence of infrared power behaved 
propagators are investigated in the Landau gauge. We also briefly develop the same ideas for the massless quark 
equation. 

THE EQUATIONS 

The QCD Dyson-Schwinger equations for the gluon and ghost propagators are (provisionally neglecting the four- 
gluon vertex and setting the number of quarks to zero): 

[IVGOP 1 = [DIM' 1 (-l)^9o 2 / -0^Gl(-r, q )A(q)G„(q-r)A(-r) (I) 

- \n c .91J ^r^(-p, 9j r)I?*ft(g)lt^(p 1 -r J -g)I? «*(r) 1 

[A(p)]- 1 = [A^p)]- 1 - N c glj -0^Gl(p,q)A(q)G v (q,p)D^(r), (2) 

where go is the bare coupling, D^ LV the gluon propagator, A the ghost propagator, rJL, the triple-gluon vertex, G M 
the ghost-gluon vertex, the zero superscript denotes bare quantities, and r = p — q. Note that the contraction over 
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color indices has already been performed, and the extra factor (—1) in front of the ghost loop is due to the Fermi 
statistics of the ghost field. 

The general expressions for the full gluon and ghost propagators in a covariant gauge £ can be written as 



A(p) 



9^ 



iG{p 2 



(3) 
(4) 



and we shall refer to the scalar functions F(p 2 ) and G(p 2 ) as the gluon and ghost form factors. 

Although we do not know the correct expressions for the full triple- gluon and gluon- ghost vertices, we can always 
formally rewrite Eqs. (|l|, §), for arbitrary vertices, as: 
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N c g 2 a 
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8tt 3 



A 2 r -K 



dy d9 y sin 2 9 [M(x, y, z) G{y)G{z) + Q(x, y, z) F(y)F(z)] 
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A 2 m 



Jo 



dy d9 y sin 2 9 T{x, y, z) G{y)F{z) . 



(5) 
(6) 



where we first contracted Eq. (W) with —iV ti,y (p)/3p 2 , where V^ljp) — g^ v — ip^p^/p 2 , thus avoiding spurious 
quadratic ultraviolet divergences iH , and multiplied Eq. (||) with i/p 2 . Subsequently we applied a Wick rotation to 
both equations, which introduces an imaginary factor i for each four dimensional momentum integration, evaluated 



two trivial angular integrations, and set x = p 2 , y = q 2 , z 
are defined by: 



The kernels M, Q, T depend on the full vertices, and 



M(p 2 ,<?,r 2 ) 
Q(p 2 ,q 2 ,r 2 ) 
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6p 2 q 2 r 2 



V^(p)Gl{-r,q)G v {q,-r), 
V^{p) TX a2 (-p, q, r) T% 2f3i (p, -r, -q) 



g ^(q)g^(r)+a 9 T Pl (q) 



'r 2 F{r 2 ) 



q 2 F{q 2 ) 



, 2 q^q 131 r^r 13-2 
1 q 2 F(q 2 ) r 2 F(r 2 ) 



T(p 2 ,q 2 



r = - 



p2 g2 r 2 



r 2 F{r 2 ) 



G °u,(P,q) G v (q,p), 



(7) 
(8) 

(9) 



where g^{q) — g^ iV ~ q^q v jq 2 . Note that the full vertices will depend on the propagator functions F and G by the 
intermediary of the three-point Dyson-Schwinger equations. 

The full, regularized, unrenormalizcd Green's functions are potentially divergent as we take the ultraviolet cutoff 
A to infinity, and we therefore introduce renormalized gluon and ghost fields, and a renormalized coupling. The 
renormalized gluon and ghost form factors Fr and G R satisfy the following multiplicative renormalization condition: 



F(p 2 ,A 2 ) = Z 3 (n 2 ,A 2 )F R (p 2 ,n 2 ), 



G(p 2 ,A 2 )=Z 3 (v 2 ,A 2 )G R (p 2 ,fi 2 ), 



(10) 



where F R (/j, 2 , /i 2 ) = Gfl(/i 2 ,/i 2 ) = 1 and Z3, Z3 are the renormalization constants for the gluon and ghost fields. 
Because of gauge invariance the renormalized coupling g satisfies 



Km 2 ) 



~, V,a 2 ) Z 3 1/2 ( m 2 ,A 2 )Z 3 (m 2 ,A 2 ) , a3 . 



(11) 



where Z\, Z\ are the renormalization constants for the triple-gluon and gluon-ghost vertices. From Eqs. ( ^o| , [ll|) it is 
easy to see[fl2| that the product 



a(x, A 2 ) 



4tt 



Zf(^,A')F R (x, ^)G R (x, fi 2 ) = &(jS,A')F R (x, ^)G r (x, tf) 



(12) 



is renormalization group invariant, i.e. independent of /i. 

We know from Taylor [fl5| t hat Zi([j, 2 , A 2 ) = 1 in the Landau gauge, and there the renormalization group invariant 
quantity a(a;,A 2 ) of Eq. (p_2[) is independent of A 2 , and nothing else but the running coupling a(x) = g 2 (x)/4:TT. 
Although all our formal derivations will be valid for any covariant gauge, in the following we will focus on the Landau 
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gauge and systematically replace a(x, A 2 ) by a(x), keeping in mind that a(x, A 2 ) is the correct renormalization group 
invariant quantity to be used when studying the equations in an arbitrary covariant gauge. 

After introducing renormalized quantities in Eq. (^) using Eqs. ([ll], [O, 12), we find the following equation for the 
renormalized ghost form factor: 



G R (x) 



z 3 (^ 2 ,a 2 



2tt 2 



dydQysm 2 9T(x,y, z) G R (y)F R (z) , 



(13) 



where a M denotes the value a(/i 2 ) of the coupling at the renormalization point. 

In an analogous way we derive an equation for the renormalized gluon form factor from Eq. 



Fr(x) 



W,A 2 



8tt 3 



dy d6 y sin 



(14) 
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Zf(n 2 , A 2 )M(x 7 y, z)G R {y)G R {z) + Z 2 ([i 2 , A 2 )Q(x, y, z)F R {y)F R {i 



From Eq. (O) we see that: 



Zi(» 2 ,A 2 ) 



Z 3 (^ 2 ,A 2 ) 



W,A 2 ) 

and after using Eq. ( |l0| ) we can write this as: 

F(v,A 2 ) G R (w,^ 



^i(/AA 2 ), 



F R (v,n*) G(w,A 2 ) 



Zi([i , A ), for arbitrary v,w < A. 



(15) 



(16) 



This expression enables us to factorize Z\ in a fi and a A dependent part in the Landau gauge [Z\ = 1). Although 
Z\ depends only on [i and A, we can always write it as the ratio ( |l6|) with arbitrary momenta v and w. We eliminate 
the factor Z\ in Eq. ( fL4| ) by substituting Eq. ( |l6| ) twice, once choosing v = w = y, and once with v = w = z. Hence, 
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Fr(x 

where we have introduced the notation 



Z^ 2 ,A 2 ) Nca » 



2tt 2 



dy dd y sin 2 9 R(x, y, z) G R {y)G R (z) 



R{x 7 y,z) = M(x,y,z) 



F(y,A 2 )F(z,A 2 ) 
G(y,A2)G(z,A*) 



Q(x,y,z). 



(17) 



(18) 



When introducing the renormalized coupling g(p 2 ) in the gluon loop of Eq. ( |l4| ) we found it most logical to use 
the Zi-identity of Eq. (p"T|). An equivalent, more direct, way is to use the Zi-identity instead. Obviously the same 
equation, Eq. jl^), results. 

The two remaining renormalization constants Z% and Z3 in Eqs. (|l^, can now easily be eliminated by subtracting 
each of the equations at two different momenta: 
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dy d8 y sin 



R(x,y,z)G R (y)G R (z) ) - ( x 
T(x,y,z)G R (y)F R (z 



x <-> a 



(19) 
(20) 



where the momentum a can, but docs not have to be taken to coincide with the renormalization point /x. Let us 
now have a closer look at Eqs. (|l^, ||(]). First of all we observe that, in the Landau gauge, the equations are free 
of renormalization constants, and that the only renormalization point dependence comes through the renormalized 
coupling and the renormalized propagator functions F R (*, /i 2 ) and G R (*, /i 2 ). Indeed, as can be seen from Eqs. (0, 
[|, ||, 18), the unspecified kernels R(x,y,z) and T(x,y,z) are to be calculated using the triple-gluon and ghost-gluon 
vertices, and the gluon and ghost form factors, all in their full, regularized, unrenormalized form, and are therefore 
independent of /i. 

On the other hand, we also know that, except for the upper integration limit, the kernels R and T are the only A- 
dependent factors remaining in the equations in the Landau gauge. If multiplicative renormalizability is to be a general 
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nonperturbative feature of gauge theories, then F R (x) and G R (x) will remain finite as A is taken to infinity. Therefore 
the kernels of the integrals in Eqs. ([19| |2(J) are such that after subtraction and integration the A-dependence vanishes 
as A is taken to infinity. Moreover, the A-dependences of the integrals in the unsubtracted equations, Eqs. ([D| p~7| ) , 
correspond to those of Z3 and Z% known from perturbation theory. This can easily be verified from the A-dependence 
of the various elements in the equations given by the renormalization group equations. Any relevant approximation 
to the propagator equations will obviously have to respect this property. 

Let us show how the multiplicative renormalizability of the solutions to Eqs. (E9[ ^Oj) is satisfied. From Eq. ( |l0| ) we 
see that for solutions to be multiplicatively renormalizable they have to satisfy 

ip 1 2n F R (x,fi 2 ) 2 G R (x,fi 2 ) 

F « {X ' V)= F R [y*,fY GR{x ' v) = G R {v>,fY (21) 

when renormalized at different scales and v 2 . It is straightforward to see, using Eq. ( p"2|) , that if F R (x,fj?) and 
G R (x, u 2 ) are solutions of the set of equations fll9| , |2(i| ) renormalized at /x 2 , then indeed F R (x~v 2 ) and G R (x, v 2 ) defined 
by Eq. (|2l| ) will satisfy the same equations now renormalized at v 2 . It is precisely using this reasoning that one can 
show that the approximations introduced in Refs. JTT| |l2| , [D| to treat the gluon loop are violating multiplicative 
renormalizability. Note that any momentum scale can be chosen as renormalization point; it is not restricted to the 
small coupling region of perturbation theory. 



Of course, the new gluon equation, Eq. (17), is completely equivalent to the original one, Eq. (|14|), both having 
the same multiplicatively renormalizable solutions. So what are the advantages of the new formulation? In contrast 
to previous studies, we have eliminated the renormalization constant Z\ without introducing any approximation 
yet, and therefore multiplicative renormalizability is preserved. We will show how this allows us to devise tractable 
approximations to the equations, having solutions that will still respect multiplicative renormalizability. With the new 
approximation we will then try to answer the following questions. Do the new nonperturbative solutions agree with 
perturbation theory for high momenta? Do the solutions behave as power laws in the infrared, and is it indeed so that 
only the ghost loop contributes to the leading infrared power in the gluon equation as deduced in Refs. |n], 12, [l3[ ? 



In fact, from Eqs. ([L9| |2fj) one can introduce novel truncations to the equations for the renormalized gluon and ghost 
propagator that will readily contradict this statement. Although we know that the bare vertex approximation from 
Ref. |L2] does not preserve multiplicative renormalizability, we can now define an alternative bare truncation using 
Eqs. ( p_9| , |20| ) that will preserve both multiplicative renormalizability and the correct perturbative limit by simply 
taking the bare approximation to the complete kernels R and T. The bare approximation is defined by substituting 

r^ 7 = rX G^ = Gl, F = F = 1, G = G = l (22) 

in Eqs. (| || |l§), hence yielding 

R(x, y, z) w Rq(x, y, z) = M (x, y, z) + Q (x, y, z) , T(x, y, z) w T (x, y, z) . (23) 



The expressions for Rq and Tq are given in Eqs. ([Al|, A2). For the gluon loop the approximation amounts to 



assume that the ratio of unrenormalized form factors F(y)F(z)/G(y)G(z) in Eq. ( |18[) cancels the full unrenormalized 
corrections to the bare triple-gluon vertex to a good approximation, at least after contraction and integration. In fact 
this could be considered as an extension to the approximation introduced by Mandelstam^]. There the Ward identity 
was used to simplify the product of full vertex and full propagators. Our study now also includes the renormalization 
aspects of the gluon loop, which were omitted before. We propose to analyze the validity of our assumption in a study 
of the 3-point Dyson- Schwinger equation in future work. 

In Section (IX) of Ref. Jl^| it was shown that consistency with perturbation theory in the high momentum region is 
not trivially satisfied in a naive treatment of the gluon loop: even if the truncation leaves invariant the perturbative 
expansion of the equations, the high momentum limit of the nonperturbative solutions usually does not agree with 
the results of perturbation theory. However, the same study pointed out that, in the bare vertex "ghost- loop only" 
truncation, the structure of the loop integrands ensuring multiplicative renormalizability is exactly what is needed 
such that the high momentum limit of the nonperturbative solutions be consistent with the results of perturbation 
theory, i.e. it yields the correct values for the anomalous dimensions of the gluon and ghost propagators, and for the 
leading order beta coefficient of the running coupling. All these arguments now also apply to the bare approximation 
to Eqs. which now include the gluon loop, and the nonperturbative solutions do agree with the leading order 

perturbative results in the high momentum region, as is illustrated in more detail in Appendix |a|. The nonperturbative 
solutions have a leading ultraviolet logarithmic behavior 

s 47T 

F R (x)~\og<x, G R (x) -log 5 x, a(x)~ — — , (24) 

f3 \ogx/A z QCD 
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with Pq = HJV c /3, 7 = — 137V c /6/3o = —13/22, and 5 = —3N c /4[3q — —9/44, thus reproducing the correct perturbative 
results for the anomalous dimensions of the gluon and ghost propagators, and for the leading order beta coefficient 
of the running coupling. Note that by satisfying the perturbative limit, the approximation also ensures that the 
unsubtracted equations Eqs. ( |l3| , |l7j ) yield the correct A-dependence for Z$ and Z%. 

Next we analyze the infrared behavior of Eqs. (^, p0[ ) in our bare approximation. Following the arguments of 
Ref. Jl3| we deduce that, as before, both gluon and ghost equations are individually satisfied by propagators which 
behave as power laws in the infrared, 

F R (x) ~ x 2k , G R (x) — x~ K , a(x) — > constant , (25) 

and, because of Eq. (^) , these power laws lead to an infrared fixed point for the running coupling, 

a = lim a^F R {x)G 2 R {x) . (26) 

x— *0 

After substitution of the power laws ( p5| ) in Eqs. (|l9| ^), the ghost and gluon equations each yield a relation 
between the infrared fixed point ao and the leading infrared exponent k: 

N c a = — ^-r , N c a = — ^ , (27) 

Xgh(n) XglW 

which are derived by equating the coefficients of the leading power of x for i->0on both sides of each equation, and 
where X 9 ft,(rt), Xgi( K ) are computed by solving the loop integrals as detailed in Ref. Jl3|]. A consistent infrared power 
solution requires the gluon and ghost equations to be satisfied simultaneously, 

Xgh{n) = Xglin) ■ (28) 

The solution of this equation yields the value of the leading infrared exponent k, and the corresponding ao can then 
be computed from Eq. (p7|). 

It is important to note that, in contrast to the studies pi [l2| the gluon loop now also contributes to the leading 
order infrared power of the gluon equation, Eq. (|l9|), and hence to Xgl{ K )i m our novel truncation. The expressions 
for Xgh{ft), Xgi( K ) resulting from the leading infrared analys is in t he b are approximation are expressed in terms of 



generalized hypergeometric functions, and are given in Eqs. ( All , A12). Numerical evaluation, using Mathematica, 
has shown that Xgi( K ) < f° r an Y K f° r which the integrals converge. This disagrees with the ghost equation 
for which Xgh{n) > 0, and hence no consistent solution for k exists, and power laws are not infrared asymptotic 
solutions of Eqs. (pjL 20) in the bare approximation. Note that in previous studies pi [l^, ^3|, in which the gluon 
loop was subleading and only the ghost loop contributed to the gluon vacuum polarization, the power law solutions 
were consistent and the value of k varied between 0.7 and 0.9 depending on the truncation. In our study the gluon 
equation seems to reject the power law solutions, and although this could at first be seen as a setback, we do not 
believe that this invalidates the general ideas of our approach. We rather believe it is an artifact of the truncation, 
and more care has to be taken in the analysis of the various contributions to the gluon vacuum polarization, as we 
will show below. A case can be made for the eventual survival of infrared power laws by observing that Xgh{i^) 1 
determined from the ghost equation, can be identified with the reciprocal of the infrared fixed point of the running 
coupling (using Eq. (|27j)), and that its numerical values can be considered physical for n in the convergence region of 
the integrals, e.g. it varies from 0.87 to 0.24 for n varying from 0.2 to 1. 

Several contributions to the gluon vacuum polarization were neglected in the gluon equation, Eq. ([!]), and we now 
extend the scope of our investigation by considering its complete form: 

[D^p)}- 1 = [DIM] 1 7$(p) - - 7$ (p) - 4l(p) - 4l d (p) - , (29) 

including the contributions from the ghost loop, gluon loop, three-gluon loop, four-gluon loop, tadpole diagram and 
quark loop. First, we will show that massless quarks can contribute to the leading infrared power analysis, and how 
they can affect the existence of a consistent infrared behavior. Although quark masses will most probably alter these 
conclusions, the main purpose of this exercise is to illustrate, in a simple way, that the absence of a consistent leading 
infrared power k, resulting from the multiplicatively renormalizable treatment of the gluon loop, does not have to 
be definitive; other contributions, overlooked till now, may reverse this situation. In fact, the hope is that this will 
be realized by the three-gluon and four-gluon loops, both involving the four-gluon vertex. Although no numerical 
calculation could yet be performed, we will schematically show that these loops can also contribute to leading order 
in the infrared, even though they are only subleading in the ultraviolet. 

Consider the quark loop contribution to the gluon vacuum polarization in Eq. (p9|): 

(30) 



KM = (-i)^ 2 / ^Tr[T'lf' (-r, q ,-p)S F ( q )r^( q ,-r,p)S F (-r)] 
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where Sf is the quark propagator, TJf is the quark-gluon vertex, the factor (—1) signifies that the extra minus sign 
is due to the fermionic character of the quark loop, and the trace over color indices has already been performed. 
The most general expression for the full quark propagator can be written as 

(31) 

where we shall refer to Z(p 2 ) as the quark form factor, and M(p 2 ) is the mass function. 
Because of gauge invariance and multiplicative renormalizability we know that 

2 zl'\p\A 2 )Z 2 {p 2 ,A 2 ) Z(p 2 ,A 2 ) = Z 2 ( f , 2 ,A 2 )Z R (p 2 , f J), (32) 

Z lf [^,A 2 ) 

where Z\ / , Z 2 are the renormalization constants for the quark-gluon vertex and the quark field, and Z R is the 
renormalized quark form factor, with Z R {p? , pr) = 1. 

The contribution of the quark- loop to Eq. ( O ) for the renormalized gluon form factor is computed from Eq. d3C|), 
after contraction, Wick rotation, and multiplication with Z%. It can always be written as 

n q R (x)=-N f ^p-Z* f (» 2 ,A 2 ) I I dyd6y S m 2 6V(x,y,z)Z R (y)Z R (z), (33) 



where we introduced renormalized quantities using Eq. (f52|), and 

V(p 2 , q 2 , r 2 ) = - ^ 2 * 2 ^ V^(p)Tr [T^°(-r, q, -p) q ■ 1 T^(q, -r,p) r ■ 7 ] , (34) 

in the massless case. All the information about the full quark-gluon vertex is contained in the kernel V . Note that 
the minus sign in front of 7r* in Eq. (^9|) has been absorbed in the definition of -k%. Similarly to the treatment of 
the gluon loop, we eliminate the renormalization constant Z\f using Eqs. (|l(], O, [32]), expressing the equality of the 
renormalized quark-gluon and ghost-gluon couplings, and find: 



Th(*) = -?5 / / dydey sin 2 6 
Z7T Jo Jo 



G R (y)G R (z) . (35) 



Z(y,A 2 )Z(z,A , s 
f G(y,A 2 )G(z,A 2 ) v ^' y ' z 



io Jo 

Hence, Eq. ([It]) remains valid, provided we replace the definition (|l^) for the kernel R by 

«* » •> - »^ -» + S£&> + 1 wMM v{ - »• " ■ 

As before the y, dependence enters only through the renormalized coupling and the renormalized gluon and ghost 
form factors, and the A dependence of V is such that the integrals are finite. We herein study the case of Nf massless 
quarks, thereby avoiding the complications introduced by the presence of explicit and dynamically generated mass 
terms. Again we introduce the bare approximation Rq to Eq. (p6|), defined by Eq. ( p2] ) supplemented with 

VI 9 = r**° , Z = Z = 1 , (37) 

such that Rq — Ma + Qo + Vq, with the quark loop component Vq given in Eq. ( A13| ). 

The ultraviolet analysis can be done completely analogously to that presented for the quarkless case in Appendix [A|, 
and it again leads to the high momentum behavior shown in Eq. (p4|), but now with (3q = (lliV c — 2iV/)/3, 7 = 
(— 13A c /2 + 2A/)/3/?o, and S = — 3iV c /4/3o, thus correctly reproducing the leading order perturbative results. 

From the structure of Eq. ( |35| ) it is clear that massless quark loops will also contribute to the leading infrared 
power in the gluon equation in our bare approximation, in the assumption of infrared power behaved gluon and ghost 
form factors. The consistency condition (|2q) is now modified as the gluon equat ion gets an extra term from Eq. (|3^), 
proportional to Nf times the quark loop contribution X q g ii K ) given in Eq. (A14): 

N cXg h(K) = N cXg i(K)+N f x q gl (K). (38) 



We analyze the effect of the quark loop contribution on the right hand side of Eq. (B8|), and observe that its 
contribution is positive. In Fig. (|l]), we show how it increases as we increase the number of massless quarks, and how 
it eventually compensates the negative value from the gluon loop (with N c = 3) such that Eq. (p8|) be satisfied. There 
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K 

FIG. 1: Determination of the infrared exponent k as the intersection of the left and right hand sides of the consistency 
condition, Eq. (^), computed from the ghost equation (dashed line) and gluon equation for Nf = 6, 8, 9, 10 (full lines). 



is a critical value Nf « 8.9 for which a consistent infrared behavior can first be found: n rs 0.76. Then, when we 
increase Nf further, two possible values of k satisfy the consistency condition, Eq. (|38|). For Nf — 9 we find re = 0.69 
and 0.82, for Nf = 10: re = 0.49 and 0.90, for Nf = 11: re = 0.37 and 0.93, for N f = 12: re = 0.28 and 0.94. As we 
keep increasing Nf — > oo, the two solutions for re go to and 1, the corresponding values for cxq are then and 47r/3. 

We conclude that, in the bare approximation, massless quarks do contribute to the leading infrared power of the 
gluon equation, and that a consistent power behavior can be found if the number of quarks is sufficiently large 
(Nf > 9), thus potentially restoring the existence of infrared power solutions for the gluon and ghost form factors. Of 
course, we know that the real world QCD does not contain such a large number of massless quarks, and the results of 
the infrared analysis will be altered, even in the bare approximation, by the presence of explicit quark masses. Further 
study is needed to investigate how massive quarks, with either explicit or dynamically generated quark masses, will 
influence these results. This will be done in future work in conjunction with a simultaneous solution of the quark 
Dyson- Schwinger equation. Note however, that the above-mentioned results might be of importance in the study of 
Yang-Mills Grand Unified Theories, which could contain a large number of massless quarks. 

More importantly, the previous calculation enabled us to see that, even though the multiplicatively renormalizable 
treatment of the gluon loop seemed to disprove the existence of infrared power law solutions, these can be reinstated 



by additional contributions to the vacuum polarization. As shown in the gluon equation, Eq. (29), even the pure gauge 
theory, without quarks, contains additional diagrams. The three-gluon and four-gluon loop contributions have been 
neglected until now, based merely on perturbative arguments. We now schematically show that this is not justified, 
and that these diagrams can also contribute to the leading infrared behavior of the gluon and ghost propagators. 
The three-gluon contribution to the unrenormalized gluon equation, Eq. (|29|), is 

41(P) = -fgoj 7^4 J T^Ce^ast-^M^ 1 ^^ (39) 

where r = p — q — k. After contraction, Wick rotation, and multiplication with Z 3 of Eq. (|39|), the additional 
three-gluon contribution to Eq. ([17]) for the renormalized gluon form factor can be written in the form: 

4 s Or) = -ZzNlgtj -0^ J ^Q 3 (p, q ,k)F(z 1 )F(z 2 )F(z 3 ), (40) 

where Zi, Z2,z 3 are the squared momenta of the gluons in the loop. Introducing the renormalized coupling using 
Eq. (pd|), we can rewrite this as 



3g/ \ _ j\r2„4^,,2^ ^1 f d q f d k 



?(x) =-AT c V(M 2 )^ J J j^Q 3 (p,q,k)F(z 1 )F(z 2 )F(z 3 ) . (41) 
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As before, we replace all the renormalization constants by the ratios of the corresponding unrenormalized and renor- 
malizcd form factors, and find 



d A q 



d A k 



(2tt) 4 J (2?r) 



F(z 2 )F(z 3 )Q 3 (p,g,k) 
G{z 2 )G(z 3 ) G 2 (z!) 



F R {z 1 )G 2 R {z 1 )G R {z 2 )G R {z 3 ) . (42) 



Using Eq. ( |l2| ) we recognize the running coupling at momentum z\, a{z\) = a ll F R (zi)G R {zi) , and can therefore 
rewrite this as 



4»=-16^ 2 iV c 2 



d*k 



(2tt) 4 J (2tt) 



^2)^3) Qs(p,q,k) 



G(z 2 )G(z 3 ) G 2 ( Zl ) 



a{z\) G R {z 2 )G R {z 3 ) . 



(43) 



One of the effects of our manipulations has been to change one of the constant renormalized couplings in front of 
the integral into a momentum dependent running coupling a(z%) inside the integral. In the assumption of infrared 
power laws, Eq. (Eq), the running coupling has an infrared fixed point, and hence, the infrared structure of the 
integrand in Eq. ([S]) is very similar to that of Eq. ([l7|), and tt r 9 (x) will potentially contribute to the leading infrared 
power analysis of the gluon equation. 

We follow the same steps in the treatment of the four-gluon loop. Its contribution to Eq. (pit) is: 



d 4 k 



(2tt) 4 J (2tt) 



(44) 



x rXa 2a3 (-P, ?, k, r)D a ^ (q)D a ^ (k)D a ^ (r)!*^ A (-q,s, -k)D a ^ {s)Y% sPi (p, -r, -s) , 



where r = p — q — k, s = q + k. After contraction, Wick rotation and multiplication with Z 3l the four-gluon contribution 
to Eq. (TPa) for the renormalized gluon form factor can be written as 



4 9 (x) = -Z 3 Nlg\ 



Q 4 (p, q, k)F( Zl )F{z 2 )F(z 3 )F(z 4 ) . 



In analogy with Eq. ( f43| ) , we find 
ir R 9 (x) = -16n 2 N^ 



d A q f d A k 
(2^7 (2^)* 



d 4 k 

F{z 2 )F{z 3 ) Q4p,q,k)F(zi) 



G{z 2 )G{z 3 ) G 2 ( Zl ) 



a(zi)G_ R (z 2 )G_ R (x 3 ) . 



(45) 



(46) 



and this diagram too can potentially contribute to the leading order infrared power. 

An infrared power law analysis of the gluon equation including the three-gluon and four-gluon contributions, 
Eqs. ([43|, Eft), will yield a condition 



(N c a ) Xg i(K) + {N c a a ) 2 X ]f 9 {K) = 1 , 



(47) 



for the pure gauge theory, where Xgf ' 49 is to be calculated from Eqs. J43|, |46| ) . This condition replaces the second 
relation in Eq. (^7j), and is now quadratic in N c ao- The solution of Eq. ( [f7| ) for N c ao, together with the first part of 
Eq. (p7j), will yield a new consistency equation 



(48) 



replacing Eq. (|2§|). The contributions of Eqs. (|43|, |46j) involve the computation of two-loop integrals, and the approx- 
imations to the kernels have to be devised and investigated carefully. Note that these kernels too are independent of 
/j, and are such that all integrals in the subtracted equations are finite. It is crucial for the survival of the power law 
solutions as infrared asymptotic behavior of the ghost and gluon propagators, that a leading infrared exponent k is 
found which satisfies the condition, Eq. (|4^). A more detailed analysis of Eqs. |^|) will be performed in future 
work. 

Finally, we briefly turn our attention to the quark Dyson-Schwinger equation, in the massless case, and examine 
how it transforms in our new approach. The equation is: 



where C F = (A 2 



[SfIp)]- 1 = {sUpT 1 

l)/2N c = 4/3, for N c = 



C F gl 



d 4 g 
(2^) 4 



17.° (p, q , r )S F (q)Tl 9 (q, p, -r)D» v (r) , 



(49) 
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After substituting Eq. (|3lj) in Eq. (J49j) , we can derive the following formula for the quark form factor: 

.2 /-A 2 F Tt 

dy d9 y sin 2 9 U(x, y, z)Z{y)F{z) , 



(x) 8tH» J Q J 



Z( 

where the kernel U depends on the full quark-gluon vertex. Multiplying Eq. 
introduce the renormalized coupling, we obtain 



(50) 



by Z 2 , and using Eq. (ttlh to 



1 



= Z 2 (» 2 ,A 2 )- 



C F g 2 (fi 2 ) Z 2 Z 2 
8tt 3 



z r( x ) °"" ^3^3 Jo Jo 

As before, we remove all the renormalization constants in the self-energy term 



dy d9 y sin 2 9 U(x, y, z) Z{y)F(z) . 



1 



Zr{x) 



Z 2 ^ 2 ,A 2 



2tt 2 



Jo 



dy d6 y sin 



z 2 (y) 

G 2 (z) 



U(x,y,z) 



Fr{z)G%{z) 



z R {y) 



Similarly to the derivation of Eq. (fl3|), we recognize the running coupling a{z) and find: 



Jo 



dy dO y sin 2 



G 2 (z) 



U(x,y,z) 



a(z) 



Z R (y) 



(51) 



(52) 



(53) 



Subtraction of this equation at two momenta allows us to eliminate the remaining renormalization constant Z 2 . Note 
that the multiplicative renormalizability of the solutions is manifested in a very simple way by the I/Zr structure 
on both sides of the equation. In fact, the shape of the equation is somewhat reminiscent of the massless, quenched 
QED, fermion equation derived by Curtis and Pennington jl6| using a vertex extension |L7| to the Ball-Chiu vertex Jl8[, 
constructed using conditions of multiplicative renormalizability. Although we leave a detailed study of the quark 
equation to future work, we already see that we can introduce a bare approximation to Eq. (|53|), in analogy to the 
procedure applied to the gluon equation, which will yield solutions respecting multiplicative renormalizability. 

Furthermore, we anticipate that, in a similar way, the scalar part Al(p 2 ) of the quark propagator will also be 
driven by the running coupling under the integral, and it will be interesting to study how dynamical chiral symmetry 
breaking gets realized with a kernel proportional with a(z), which has an infrared fixed point if the infrared power 
laws, Eq. (p5|), are verified. This will also be investigated in more detail in a further publication. 

Note that we have used the running coupling a(x) = a^F^x, fi 2 )Gj t (x, fi 2 ) in the derivations of Eqs. ( |43] , 
fl6| , |53| ), but we could as well introduce the more general renormalization group invariant quantity a(x,y,z) — 
&fj,Fii(x, /i 2 )Gfi(y, /j, 2 )Gji(z, /i 2 ). Although this will change the appearance of the equations, it will not change their 
solutions, as long as no approximations to the kernels are introduced. However, when devising truncations to the 
kernels, the choice of momenta should be handled judiciously for the approximations to make sense. 



CONCLUSIONS 



We have reformulated the coupled set of continuum equations for the renormalized gluon and ghost form factors in 
QCD, such that the multiplicative renormalizability of the solutions is manifest, independently of the specific form of 
full vertices and renormalization constants. In the Landau gauge, all the renormalization constants are eliminated, and 
the renormalization point dependence is completely transposed into the renormalized coupling and the renormalized 
form factors. 

Furthermore, the specific structure of the equations allows us to devise novel approximations, necessary to make 
them tractable, with solutions which respect the principles of multiplicative renormalizability and have the correct 
leading order perturbative limit. This had not been achieved in a consistent way before. 

We showed that the gluon and ghost equations are each individually satisfied by infrared power behaved gluon and 
ghost propagators. However, these infrared power solutions can only truly be identified with the QCD propagators 
if both equations are satisfied simultaneously, thereby determining the value of the leading infrared exponent and of 
the infrared fixed point of the coupling. 

In contrast with previous research, our study shows that the leading infrared behavior of the gluon equation is 
not solely determined by the ghost loop. The new approach shows that the contributions of the gluon loop, the 
three-gluon loop, the four-gluon loop and even of massless quarks should be taken into account when performing the 
leading order infrared analysis. 

Moreover, we showed that, in our new Landau gauge truncation, the gluon loop contribution to the gluon vacuum 
polarization removes the existence of a consistent infrared power solution. It is our hope that additional contributions 
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to the gluon vacuum polarization will reinstate the power law solutions. As an illustration we showed how including 
a large number of massless quarks is one way to recover consistency. The three-gluon and four-gluon loops will 
also contribute to the leading infrared power behavior in this approximation, but its quantitative treatment is more 
complicated and requires further investigation. We are therefore not yet able to conclude whether or not the gluon 
and ghost propagators are power behaved in the infrared. 

We have also briefly shown, for the massless case, how the quark equation can be treated in a similar way, and 
how the specific structure of the equation will allow us to construct truncations with multiplicatively renormalizable 
solutions for the quark propagator. 
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APPENDIX A: BARE APPROXIMATION 



The bare approximation to Eqs. (M, EG), defined by Eqs. yields the following kernels: 



R (x,y,z) 



8y 2 z 2 

-(-- 

ix \y 2 

--2 + 
V 



+ JL /1 

v z \y 
n n 

y z 

y 

X 









) Az 2 + 


G- 





1 /15 
8 {tf 
1 

' 2X 2 
1 

2~ 



34 

yz 



15 



z 2 6z 

H 14+ — + 



sin 



Axy 



6y 

z 



(Al) 



(A2) 



Ultraviolet behavior 

We follow the arguments of Ref. [jl2| Section (VII)] to study the ultraviolet behavior of Eqs. ([H], 20) with kernels, 
Eqs. (Al, |A2| ). The leading ultraviolet order solutions satisfy 

2 

T^T = f dy-Gl(y), (A3) 

Fr{x) 4tt J x y 

TT^T = 1"? / dy^F R (y)G R (y), (A4) 
G R {x) 4tt J x y 

where fo and go are easily computed from the kernels. We have chosen the renormalization scale \x as subtraction 
point, and both x and fj, are in the perturbative region. It is easy to see that this is solved by: 



F, 



where (3q = fo + 2go, 7 = —fo/Po, and 5 = —go/Pa- This leads to a running coupling: 



a(x) = a Mx) GUx) = £ = ^^r~ • 

47T & Tt 1 



A 5 " 



where we define 

A? = ii 2 pm I - 



a2 qcd = ^ 2 exp ( -—J- ) (A7) 



to leading order. Hence, the solutions flA5|) can be rewritten as 



Fje(aj,/i a )=( — ) ~log 7 x, G^x, ^) = I ^ ) "log's, (A8) 
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which is in agreement with the renormalization group equation results. In the bare approximation we find (3q = 1 lN c /3, 
7 = -13iV c /6/3 = -13/22, and S = -3A c /4/3 = -9/44. 



Infrared behavior 



In this section we give the results for Xgh( K ) ancl Xgi( K ) °f Eq. (|27}). They are determined from the coefficients 
of the leading power of x for x — ► 0, after substitution of the power laws ( p5| ) in Eqs. (JH| ^0|), and solution of the 
integrals. 



The angular and radial integrals are readily solved using the method of Ref. 13 , and the ^-dependence of a typical 
integral is given by: 



/ dyy b d6 sin 2j ' z a ~ H (a, b) x a+b , 
Jo Jo 



with z — x + y — 2 Jxy cos 9, and 



H 3 {a,b) = B [j 



1 1 

2' 2 



6 + 1 



iF 2 (-a,-a-j,b+l;j + l,b + 2;l) 



(A9) 



(A10) 



1 



a + b+1 



}F 2 (-a, -a- j, -a - 6 - 1; j + 1, -a - 6; 1) 



where £? is the Beta-function, 3_F2 is a generalized hypergeometric function. The integrations are performed using the 
integration formula (3.665.2) of Ref. |19| and the definition of hypergeometric functions pG] . 



Using Eq. (A£) it is straightforward to see that the ghost equation, Eq. (EOT) , with bare kernel (A2), yields: 



Xgh(K) = -^H 2 (2k- 2,1 -k). 



(All) 



The gluon equation, Eq. (|l^), with kernel (Al), gives a similar expression 



48tt 2 



#i(-2 - «, -1 - k) - 15 Hi(-2 - At, 1 - At) + 2ifi(-2 - At, 2 - At) (A12) 

+4i?i(-2 - /t, 3 - At) + 8Hi(-2 - k, -k) + 8F X (-1 - At, -1 - n) - 22ff x (-l - At, 1 - At) 
+24 H x (-1 - At, 2 - At) - 34 H x (-1 - At, -At) + 2 - At, -1 - At) + 24 ffx(l - At, -At) 

+4if x (2- At,-1 - At) - 15ifi(-At,-l - At) - 56iii(-At, 1 - At) - 22ifi(-At,-At; 



Including massless quarks 

The bare, massless quark approximation Vq, to the kernel V of Eq. (|35j), is given by: 



Vo(x,y,z) 



2 


' 1 




(r 


4) 


-Uf- 


3 


yz 


X 









(A13) 



After substituting the power solutions, Eq. (|25j), in the quark loop, Eq. (35), and solving the angular and radial 
integrals, we find a quark contribution: 

x i gl ( K ) = — h 1 (-1-k,,1-k)-2H 1 (-1-k,2-k)+H 1 (-1-k,-k,)-2H 1 (1-k,-k) (A14) 

+4iJl(-At, 1 - At) + i/l(-Ai, -At) 



per massless quark flavor. 
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